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Abstract 

We examine an algorithm for the visualization of domain-coloured Riemann 
surfaces of plane algebraic curves. The approach faithfully reproduces the 
topology and the holomorphic structure of the Riemann surface. We discuss 
how the algorithm can be implemented efficiently in OpenGL with geometry 
shaders, and (less efficiently) even in WebGL with multiple render targets 
and floating point textures. While the generation of the surface takes 
noticeable time in both implementations, the visualization of a cached 
Riemann surface mesh is possible with interactive performance. This 
allows us to visually explore otherwise almost unimaginable mathematical 
objects. As examples, we look at the complex square root and the folium 
of Descartes. For the folium of Descartes, the visualization reveals features 
of the algebraic curve that are not obvious from its equation. 


1 Introduction 

1.1 Mathematical background 

The following basic example illustrates what we would like to visualize. 
Example 1.1. Let y be the square root of x, 

y = y/x. 

If cc is a non-negative real number, we typically define y as the non-negative real 
number whose square equals a;, i.e. we always choose the non-negative solution 
of the equation 

2/^ — X = 0 

as 2 / = y/x. For negative real numbers x, no real number y solves 
However, if we define the imaginary unit i as a number with the property that 
i^ = — 1 then the square root of x becomes the purely imaginary number 

y = i^/\x\. 
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Figure 1.2: When a complex number (black points) runs along a circle centred 
at the origin of the complex plane, its square roots (white points) move at 
half the angular velocity (left image). After a full turn of x, the square roots 
have interchanged positions. The real part of the square root that initially had 
positive real part has become negative, and vice versa (right image). 


Together, these two conventions yield a continuous square root function 

yi: C. 


For complex numbers x, Equation 1 has exactly two complex solutions (counted 
with multiplicity), the square roots of x. 

We have seen that, for real numbers x, we can choose one solution of |Equation 1| 
for the square root and obtain a square root function that is continuous over the 
real numbers. In contrast, we cannot for every complex number x choose one 
solution of [Equation l| so that we obtain a square root function that is continuous 
over the complex numbers: If we plot the two solutions of [Equation Ij as x runs 
along a circle centred at the origin of the complex plane, we observe that y moves 


at half the angular velocity of x (see Figure 1.2). When x completes one full 
circle and reaches its initial position again, the square roots have interchanged 
signs. Therefore, a discontinuity occurs when x returns to its initial position 
after one full turn and the square root jumps back to its initial position. 

Note that by choosing the values of the square root in a different manner, 
or, equivalently, letting x start at at a different position, we can move the 
discontinuity to an arbitrary position on the circle. Moreover, note that there is 
(at least) one discontinuity on any circle of any radius centred at the origin. 

In order to define the principal branch of the complex square root function, 
we usually align the discontinuities along the negative real axis, the canonical 
branch cut of the complex square root function, and choose those values that on 
the real axis agree with the square root over the real numbers. 

Alternatively, we can extend the domain of the complex square root to make 
it a single-valued and continuous function. To that end, we take two copies of 
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the extended complex plane and slit them along the negative real axis. On the 
first copy, we choose the solution of Equation 1| with non-negative real part as 
the complex square root of x; on the second copy, we choose the other solution. 
We glue the upper side (lower side) of the slit of the first copy to the lower 
side (upper side) of the slit of the second copy and obtain a Riemann surface 
of the complex square root. (In three dimensions, this is not possible without 
self-intersections.) 


Remark 1.3. On the Riemann surface, the complex square root is single-valued 
and continuous. It is even analytic except at the origin and at infinity, which 
are exactly the points where the two solutions of Equation 1 coincide. 


Remark 1.4. The branch cut is not a special curve of the Riemann surface. When 
we glue the Riemann surface together, the branch cut becomes a curve like every 
other curve on the Riemann surface. If we had used a different curve between 
the origin and infinity as branch cut, we would have obtained the same result. 


Remark 1.5. [Equation 1| describes a parabola. We can proceed analogously to 
obtain Riemann surfaces for other plane algebraic curves. 


1.2 Previous work 

Probably the most common approach for visualization of functions is to plot a 
function graph. However, for a complex function 

5: C^C, 


the function graph 


{{z,g{z)) l^eCjcCxC-M^x 


is a (real) two-dimensional surface in (real) four-dimensional space. 

One way to visualize a four-dimensional object is to plot several two- or 
three-dimensional slices. This approach seems less useful for understanding the 
overall structure of the object. 

Another traditional method to visualize complex functions is domain colour¬ 
ing. The principle of domain colouring is to colour every point in the domain of 
a function with the colour of its function value in a reference image. If we choose 
the reference image wisely, a lot of information about the complex function 
can be read off from the resulting two-dimensional image (see e.g. ( [Poelke and 
Polthier, 2012) and (Wegert, 2012)). The idea of lifting domain colouring to 


Riemann surfaces is due to Poelke and Polthier (2009). 


We can interpret a Riemann surface of a plane algebraic curve 


C: f{x,y) = 0 


( 2 ) 


as a function graph of a multivalued complex function, which maps every x 
to multiple values of y. If f{x,y) is a polynomial of degree n in y, there are 
exactly n values of y for every value of x that satisfy /(x, y) = 0 (counted with 
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multiplicity). Every such pair {x^y) corresponds to a point on the Riemann 
surface. In other words, the Riemann surface is an n-fold cover of the complex 
plane. 

Let TT : {x,y) 1-^ X denote a projection function on the Riemann surface. Then 
the values of y at a: correspond to the elements of the fibre 7r“^(a:). The situation 
is analogous to function graphs of single-valued functions from the real numbers 
(or the real plane) to the real numbers, where one function value lies above every 
point in the domain. 

We can transfer the Riemann surface from (real) four-dimensional space into 
(real) three-dimensional space by introducing a height function iL: C —>■ K. We 
typically use the real part as a height function. We plot the surface 

{(Rea;,Imx,iL(y)) | x,y G C: /(x,y) = 0} 


and use domain colouring to represent the value of y at every point of the surface. 

In practice, we want to generate a triangle mesh that approximates the 
Riemann surface as the graph of a multivalued function over a triangulated 
domain in the complex plane. The Riemann surface mesh approximates the 
continuous Riemann surface in the following sense: The y-values at the vertices 
of a triangle of the Riemann surface mesh result from each other under analytic 
continuation along the edges of the underlying triangle in the triangulated 
domain. If f{x, y) = 0 is a polynomial of degree n in y there are n values of y 
above every vertex x of the triangulated domain. Hence, we have to determine 
which of the 3n values of y above a triangle in the triangulated domain form 
triangles of the Riemann surface mesh. A wrong combination of values of y to 
triangles might for example occur due to discontinuity if we used the principal 
branch of the square root function for the computation of y. This would produce 
artefacts in the visualization for which there is no mathematical justification. 

For the generation of such a Riemann surface mesh, previous algorithms 
have solved systems of differential equations (Trott, 2008 Nieser et ah, 20101 or 


explicitly identified and analyzed branch cuts to remove discontinuities (Kranich, 
2012 Wegert, 2012[ Section 7.6). 


In the next section, we discuss an algorithm based on a different idea: We 
can exploit that y(x) is continuous almost everywhere on the Riemann surface 
and therefore, if x changes little, so does y(x). 


2 Algorithms 

In this section, we describe algorithms for generating and visualizing domain- 
coloured Riemann surface meshes of plane algebraic curves. Let 

C:/(x,y)=0 (3) 

be a complex plane algebraic curve. In particular, let / be a polynomial with 
complex coefficients of degree n in y. Moreover let [/ C C be a triangulated 
domain in the complex plane. (In practice, U is typically rectangular.) 
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We want to generate a Riemann surface mesh of C. The mesh discretizes a 
part of a (real) two-dimensional surface in (real) four-dimensional space. We 
can visualize it using a height function and domain colouring, as described in 
the previous section. 

We obtain a Riemann surface mesh of C as a graph of the multivalued function 
induced by [Equation 3[ which maps every value of x in [/ to n values of y such 
that /(x, y) = 0 . 

For every triangle in 17, we thus obtain n values of y at each of its three 
vertices. The problem is to determine whether, and if so, how, the 3n values of 
y can be combined to form triangles of the Riemann surface mesh. The resulting 
triangles should be consistent with the fact that y as a function of x is analytic 
almost everywhere on the Riemann surface. This is impossible if the triangle in 
U contains a ramification point of y(x). In this case, we subdivide the triangle 
to obtain smaller triangles mostly free of ramification points. Otherwise, the 
triangles of the Riemann surface mesh are uniquely determined by analytic 
continuation of y(x) along the edges of the triangle in U. 

In order to find these triangles of the Riemann surface mesh, we use the 
following idea: Consider a triangle AX 1 X 2 X 3 in U that is free of ramification 
points of y{x). Under this assumption, y{x) is continuous on those parts of the 
Riemann surface that lie above AX 1 X 2 X 3 . Hence, for every e > 0 there exists 
5 > 0 such that |y(xi) — y(x 2 )| < £ for all Xi,X 2 with |xi — X 2 I < i5. If £ is 
half the minimum distance between the n values of y(x) at xi and |xi — X 2 I is 
smaller than the corresponding S, then the values of y(x) at X 2 are closer to the 
corresponding values of y(x) at xi than to any other value of y(x) at xi. 

In other words, if triangle AX 1 X 2 X 3 is small enough, we can combine the 
values of y at its vertices to triangles of the Riemann surface mesh based on 
proximity: Among the 3n values of y at the vertices of triangle AX 1 X 2 X 3 , every 
three values of y closest to each other form a triangle of the Riemann surface 
mesh. 

We can algorithmically compute a (5 > 0 as above using the epsilon-delta 
bound for plane algebraic curves of jTheorem 3.1| [Theorem 3.1| is of essential 
importance for our approach. Our approach only works because [Theorem 3.1[ 
provides us with a reliable bound computable as a function of x that depends 
only on a few constants derived from the coefficients of /(x,y). 

If triangle AX 1 X 2 X 3 is not small enough to correctly combine the values of y 
at its vertices based on proximity, we subdivide the triangle. 

In summary, we obtain the following algorithm: 

Algorithm 2.1 (Generation of a Riemann surface mesh). Let [/ C C be a 
triangulated domain in the complex plane. Let 

C-. f{x,y) = 0 

be a complex plane algebraic curve and f{x,y) a polynomial of degree n in y. 
We prescribe a maximal subdivision depth (as a maximal number of iterations 
or as a minimal edge length). 
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1 . 

2 . 


Compute the global ingredients of the epsilon-delta bound of Theorem 3.1 
for y{x). 

For every triangle Ax 1 X 2 X 3 in U: 

(a) Compute the 3n values of y{x) at xi,X 2 ,X 3 , 

{yk{xj) I f{xj,yk{xj)) =0, j = 1,2,3, k = l,2,...,n}. 


(b) 


(c) 


Compute half the minimum distance between the values of y{x) at 
each of the vertices of Aa;ia:: 2 a; 3 , 


e{xj) = ^mm\yk{xj) - yiixj)\, j = 1,2,3. 


k^l 


Compute 6{xj) by the epsilon-delta bound of Theorem 3.1 so that 


\y{xj) - y{x)\ < e{xj), ii \xj - x\ < S{xj), j = 1,2,3. 


(d) Determine which of the edges of Aa::ia; 2 a ;3 are longer than the min¬ 
imum of the S{xj) at their endpoints and must be subdivided. 

(e) Select the right adaptive refinement pattern (see Figure 2.2) and 
subdivide Aa::ia; 2 a ;3 accordingly. 


3. Repeat step 2 until the maximal subdivision depth is reached. 

4. Discard every triangle in U with an edge longer than the minimum of the 
S{xj) at its endpoints. 

5. For every triangle Aa;ia: 2 a ;3 in U, combine the values of y{x) at its vertices 
to triangles of the Riemann surface mesh based on proximity. More formally, 
the triangles added to the Riemann surface mesh comprise the vertices 


{xi,yk{xi)), 

(x2,argmin|yfe{a;i) - yi{x2)\), 
yi{x2) 

(x3,argmin|yfc(a;i) - yi{x3)\), 


for k = 1,2,... ,n. 

6 . Output the Riemann surface mesh and stop. 
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Remark 2.3. By construction, [Algorithm 2.1 generates a Riemann surface mesh 
that is consistent with the analytic structure of the Riemann surface of C. 

Remark 2.4. The adaptive refinement patterns used for the subdivision of 
triangles, whose edges are too long, produce a watertight subdivision. 

Remark 2.5. Step 4 of [Algorithm 2.1 produces holes around the ramification 
points of y{x). We can make these holes very small if we choose the maximal 
subdivision depth appropriately. 

For the visualization of a Riemann surface mesh, we use the following algorithm: 


Algorithm 2.6 (Visualization of a Riemann surface mesh). Let a Riemann 
surface mesh and a domain colouring reference image be given. We choose a 
height function iL: C —>■ K to transform a point on the Riemann surface mesh 
from (real) four-dimensional space to a point in (real) three-dimensional space, 

{x,y(x)) I—>■ (Rea;,Imx,iL(y(a;))). 

1. Draw the mesh that results from transforming every vertex of the Riemann 
surface mesh as above. 


2. Interpolate the value of y{x) on the transformed mesh. 

3. Assign to every point on the transformed mesh the colour in the reference 
image of the value that y{x) attains at that point on the transformed mesh. 

Remark 2.7. If we choose the real (or imaginary) part of y{x) as a height function, 
the transformation from (real) four-dimensional to (real) three-dimensional space 
becomes a projection. 

Remark 2.8. Using the real part of y{x) as a height function has the advantage 
that the visualization then contains the image of C interpreted as a real plane 
algebraic curve. It is the intersection of the visualization of the Riemann surface 
mesh in (real) three-dimensional space with the Re x-Re y-plane (the xz-plane, if 
we label the coordinate axes of real three-dimensional space such that the x-axis 
points to the right and the z-axis points upwards). 

Remark 2.9. The computation of the Riemann surface mesh by [Algorithm 2.1 [ is 
independent of the choice of height function used for its visualization. 


3 Implementation 


In this section, we discuss how [Algorithm 2.1 and [Algorithm 2.6 can be imple¬ 
mented using OpenGL and WebGL. Since WebGL targets a much wider range 
of devices, its API is more limited than that of OpenGL. Gonsequently, our 
implementation using WebGL differs substantially from our implementation 
using OpenGL. Before we discuss each setup separately, let us talk about what 
they have in common. 

The main part of our programs is written in shading language (GLSL for 
OpenGL and ESSL for WebGL) and runs on the GPU. We use the CPU to 
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compute the global ingredients for the epsilon-delta bound, to generate shading- 
language code that computes the epsilon-delta bound for C: f{x,y) = 0 as a 
function of a;, and to generate a coarse triangulation of the input domain. 

The implementations in OpenGL and WebGL share some shading language 
code. Since there is no native type for complex numbers, we represent them 
using two-dimensional floating point vectors. Common routines include complex 
arithmetic, numerical root-finding algorithms, the computation of the epsilon- 
delta bound, and domain colouring. 

The implementation of complex arithmetic is straightforward and we shall 
not go into detail about it. 

We need numerical root-finding algorithms to approximate roots of polyno¬ 
mials in order to compute values of y{x) (and to compute the global ingredients 
of the epsilon-delta bound). For instance, Laguerre’s method (Press et ah, 


2007 Section 9.5.1) and deflation (Press et ah, 2007 Section 9.5.3) or Weierstrafi- 
Durand-Kerner method ( Weierstrafi, 189T| ' Durand, 1960: Kerner, 1966) are 
well-suited. The latter may be a little easier to implement in shading language 
(due to the absence of variable-length arrays). 

For the computation of the epsilon-delta bound, we use the following theorem: 


Theorem 3.1 (cf. ( |Kranich, 2015 Theorem 2.1)). LetC: f{x,y) = 0 be 
complex plane algebraic curve, where 


n 

f{x,y) = J2^k{x)y^~'' 

k^O 

is a polynomial of degree n in y whose coefficients ak{x) are polynomials in x of 
the form 

rrik 

ak{x) = ^akix^^~K 

Let xi G C be a point in the complex plane at which neither the leading coefficient 
ao{x) nor the discriminant of f(x,y) w.r.t. y vanish. Then for every e > 0, we 
can algorithmically compute (5 > 0 such that 

\Vj{xi) - yj{x2)\ <e 

for all holomorphic functions yj{x), j = 1,2,... ,n, which satisfy f{x, yj{x)) = 0 
in a neighbourhood of xi and for all X 2 with \xi — X 2 \ < S. 

We obtain 
















where 


p < min{|xi - x \: ao{x) ■ Ay{f{x,y)){x) = 0}, 


Y := max 

3 


f^{xi,yj{xi)) 


M := 2 max ( — 
k \ao 


ao 


fyixi,yj{xi)) 
mo rrik 

l^ool - ^i| - p) > 0 , a/e := ^ \aki\{\xi\ + |p|)' 


-I 




/=0 


Ay{f(x,y)){x) denotes the y-discriminant of f{x,y), and xi, I = 1,2,..., mg, 
are the zeros ofao{x). 

Note that the computation is parallelizable since the epsilon-delta bound can 
be implemented as a function of x that depends on only a few constants derived 
from the coefficients of f{x,y). 

Instead of computing texture coordinates, which would depend on the range 
of y{x) on the input domain, we generate the domain colouring procedurally 
on-the-fly. To that end, we use a variation of the enhanced phase portrait colour 


scheme of (Wegert, 2012 Section 2.5). The reference image is shown in Figure 4.1 


We discuss the colour scheme in ISubsection 4.11 

The main difference between the implementations in OpenGL and WebGL 
is how the common routines can be combined to realize [Algorithm 2.1| and |A1 
jgorithm 2.6| 

3.1 An implementation in OpenGL 


3.1.1 Implementation of [Algorithm 2TT] in OpenGL 

Our implementation of [Algorithm 2.1 in OpenGL comprises three GLSL pro¬ 
grams, for initialization, subdivision, and assembly of the Riemann surface mesh. 
We cache the output of each program using transform feedback and feed it back 
to the next program. 

The initialization program consists only of a vertex shader, which operates on 
the vertices of the triangulated input domain. For every vertex x, we compute 
2 /fc(x), fc = 1, 2,... ,n, and d{x). 

After initialization, we run the subdivision program. The program consists 
of a pass-through vertex shader and a geometry shader. The geometry shader 
operates on the triangles of the triangulated input domain or of its last sub¬ 
division, respectively. We have access to the values of xj, S(xj), and yk{xj), 
fc = 1, 2,..., n, j = 1, 2, 3, at the vertices of each triangle AxiX 2 X 3 . We determ¬ 
ine which edges of triangle /\xiX 2 X'i are longer than the minimum of the 5{xj) at 
their endpoints. In order to subdivide these edges, we compute their midpoints 
X, and 5{x) and ykix), k = 1, 2,..., n, at the midpoints. We use the appropriate 
adaptive refine pattern of [Figure 2.2[ and output between one and four triangles 
for every input triangle. In doing so, we reuse previously computed values rather 
than recomputing them. We run the subdivision program iteratively until we 
reach the prescribed maximal subdivision depth. 


9 




















The assembly program consists of a pass-through vertex shader and a geo¬ 
metry shader. The geometry shader operates on the triangles of the adaptively 
subdivided input domain. We again have access to the values of Xj, S{xj), and 
yk{xj), fc = 1, 2,..., n, j = 1, 2, 3, at the vertices of each triangle AX 1 X 2 X 3 . For 
every triangle AX 1 X 2 X 3 , we test whether one of its edges is longer than the 
minimum of the S{xj) at its endpoints. In this case, we discard the triangle. 
Otherwise, we determine the triangles of the Riemann surface mesh by proximity 
Algorithm 2.1[ step 5) and output these n triangles. 


(see 


We also cache the assembled Riemann surface mesh using transform feedback 
so that we can pass it as input to our implementation of the visualization 
algorithm (Algorithm 2.61. 


3.1.2 Implementation of [Algorithm 2 T 6 in OpenGL 


Our implementation of [Algorithm 2.6] in OpenGL consists of one GLSL program 
with a vertex and a fragment shader. 

The vertex shader operates on the vertices of a Riemann surface mesh 
generated by our implementation of [Algorithm 2.1] We apply height function 
iJ: C —>■ K to map each (real) four-dimensional vertex 


(Rea;, Im x. Re y (x), Im y (x)) 


to a (real) three-dimensional vertex 

(Rex, Imx, H{y{x))). 


We homogenize the coordinates of this vertex and transform them using the 
model-view-projection matrix. We pass y(x) as a varying variable to the fragment 
shader. 

The fragment shader operates on the interpolated value of y(x) at a fragment 
of a pixel of the output device. We compute the colour of y(x) according to our 
domain colouring reference image. 


3.1.3 Remarks 


Using our implementation, the generation of a Riemann surface mesh takes 
little but noticeable time. The bottlenecks of the implementation are numerical 
root-finding and iterative subdivision. However, if we use transform feedback 
to cache the Riemann surface mesh and pass it to the implementation of the 
visualization algorithm, we obtain interactive performance. 

Another advantage of using transform feedback to cache the Riemann sur¬ 
face mesh is that we can easily export the data. If we additionally compute 
texture coordinates and a high-resolution reference image, we can even print our 


visualization using a full colour 3D printer (see Figure 3.2). 
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Figure 3.2: SD-printed models of domain-coloured Riemann surfaces of square 
root, folium of Descartes, and unit circle (clockwise). The merchandise coffee mug 
of DFG Collaborative Research Center TRR 109, “Discretization in Geometry 
and Dynamics”, is included in the picture as an indicator for the size of the 
models. 


3.2 An implementation in WebGL 

In order to support a wider range of devices, the WebGL API is much more 
limited than the OpenGL API. Particularly, in WebGL, geometry shaders and 
transform feedback are currently unavailable. (The WebGL 2 draft includes 
transform feedback and compute shaders.) Therefore our implementation in 
WebGL differs substantially from our implementation in OpenGL. 


3.2.1 How to replace transform feedback 


Instead of transform feedback, our implementation in WebGL uses floating 
point textures (specified in the OES_texture_float extension) and multiple 
render targets (specified in the WEBGL_draw_buffers extension). I do not claim 
originality of this approach. It is commonly used for running simulations on the 


GPU. The original idea may be due to (Crane et ah, 20071. We number the 


vertices of every mesh consecutively and pass this number (index) to the vertex 
shaders along with the other attributes. In particular, vertices that are shared 
among several triangles must be duplicated and numbered separately. Hence, 
we assume that every triangle appears as three consecutive vertices in array 
buffer storage (triangle soup). We use floating point textures essentially as we 
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would arrays of floats, indexed by vertex number. We store values corresponding 
to the fc-th vertex in the fc-th pixel of a texture. We can store up to four 
floats per pixel of a floating point texture, namely one float each in the red, 
green, blue, and alpha channel. If we need to store more than four floats, we 
use multiple render targets which allows us to colour the same pixel of several 
textures simultaneously. 

We want to store values we compute for a vertex in textures (‘transform’ 
in transform feedback). To that end, we bind the array buffers and draw the 
contents as points (as opposed to triangles). In the vertex shader, we compute 
the positions of the point with index k (in normalized device coordinates) so 
that it is rasterized as the fc-th pixel of the render target textures. Recall that 
normalized device coordinates range in [—1,1]^. For example, if h and w denote 
the height and width of the textures in pixels, we assign to the point with index 
k the position 


/ 2 ■ {k mod lu) + 1 2 • [k/w\ + I 

w h 

In the fragment shader, we compute the values to be stored and assign them as 
output colours in a specific order (as we later want to retrieve them). 

We want to read stored values for a vertex from textures (‘feedback’ in 
transform feedback). To that end, we bind the textures and an array buffer 
containing a range of vertex numbers (indices). We draw the contents of the 
array buffer as points or triangles (depending on whether we want to send the 
output to different textures or to the screen). In the vertex shader, we compute 
texture coordinates for the point with index k which allow us to lookup the fc-th 
pixel from the textures. Recall that texture coordinates range in [0,1]^. For a 
texture of height h and width w, we compute the texture coordinates 

{k mod w) + 0.5 (fc mod w) + 0.5 \ 
w ’ h ) ' 

Adding 0.5 in the numerators accounts for the fact that we want to obtain 
coordinates for the centre of a pixel in order to avoid interpolation with adjacent 
pixels. We pass the texture coordinates to the fragment shader, where we can 
use them to perform a texture lookup. 

In order to access data of a whole triangle (as in geometry shaders), we can, 
in the vertex shader, determine the indices of the other vertices of the triangle. 
For example, the point with index k is part of the triangle whose vertices have 
indices 



k — {k mod 3), k — {k mod 3) + 1, and k — {k mod 3) + 2. 

We compute normalized device coordinates or texture coordinates for all three 
indices and pass them to the fragment shader, together with the index of the 
triangle vertex currently under consideration. 
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3.2.2 How to replace geometry shaders 


We replace the geometry shader of the subdivision program of our implementation 
in OpenGL using a variation of a method proposed by |Boubekeur and Schlick] 
|(2008| ). The method works as follows: We precompute all adaptive refinement 
patterns up to a certain subdivision depth, in our case eight adaptive refinement 
patterns up to depth one (see Figure 2.2). We use barycentric coordinates to 
store the positions of the triangle vertices of each refinement pattern in an array 
buffer. Using array buffers of different lengths allows us to achieve variable-length 
output, as with geometry shaders. For every triangle of a coarse input mesh, 
we draw the triangles in the array buffer of the appropriate adaptive refinement 
pattern. We use the vertex positions of the input triangle (read from a texture 
or from uniform variables) and the barycentric coordinates of the triangles of 
the adaptive refinement pattern to compute the vertex positions of the output 
triangle. 

We can combine this method with floating point texture and multiple render 
targets as outlined above, if we number the vertices of each adaptive refinement 
pattern consecutively and store those indices together with the barycentric 
coordinates. We pass an offset as a uniform variable to the vertex shader that 
needs to be added to the indices. We draw the adaptive refinement pattern 
and increment the offset by the number of vertices in the adaptive refinement 
pattern. 

The geometry shader of the assembly program has fixed-length output. It 
generates exactly n triangles of the Riemann surface mesh per triangle of the 
(subdivided) input mesh. We can replace it with n invocations of a vertex shader, 
one for every sheet of the Riemann surface mesh. We pass the number of the 
current sheet to the vertex shader as a uniform variable. 


3.2.3 Remarks 

We cannot expect our WebGL implementation to reach the same performance 
as our OpenGL implementation. In the subdivision program, since we draw a 
different adaptive refinement pattern for every triangle of the input mesh, we 
lose parallelism. Consequently, subdivision in WebGL is much slower than its 
OpenGL counterpart. However, if we cache the assembled Riemann surface mesh 
(in textures) and pass it to our implementation of the visualization algorithm, 
we can still achieve interactive performance. 


4 Examples 

In this section, we discuss domain-coloured Riemann surface meshes for the 
complex square root function and for the folium of Descartes. Before that, let 
us explain our domain colouring reference image so that we can interpret the 
domain-coloured Riemann surface meshes. 
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4.1 Domain colouring reference image 

Recall that the basic idea of domain colouring is the following: If we want to 
visualize a complex function 


f: KcC^C, 


we face the problem that its graph is real four-dimensional. However, we can 
visualize the behaviour of the function by colouring every point in its domain 
with the colour of the function value at that point in a reference image. The 
reference image is the domain colouring of the complex identity function. 

Depending on what reference image we choose, we can read off various 
properties of a function from its domain colouring. For an overview of different 


colour schemes, we refer to (Poelke and Polthier, 2012 Wegert, 2012). 


As our reference image, we use a variation of the enhanced phase portrait 

Section 2.5). The reference image is best 


colour scheme of (Wegert, 2012 


described using polar coordinates 

re“^ = r (cos ip + i sin cp) 

of a complex number with modulus r > 0 and phase ip G [0, 27r). 

Firstly, we encode the phase at any point in the domain as the hue of its 
colour (in HSI colour space). In a square with side length 10 centred at the origin, 
we thus obtain the colour wheel shown in [Figure 4.1a| As the phase changes 
from 0 to 27r, we obtain every colour of the rainbow. Positive real numbers, 
which have phase 0, are coloured in pure red. Negative real numbers, which 
have phase tt, are coloured in cyan. Purely imaginary numbers do not have 
such distinctive colours. (This can be fixed using the NIST continuous phase 
mapping, which scales the phase piecewise linearly so that purely imaginary 
numbers with positive imaginary part become yellow and purely imaginary 
numbers with negative imaginary part become blue. See (NIST, 2014 http: 
//dlmf.nist.gov/help/vrml/aboutcolor#S2.SS2) 
follow this approach here.) 


For simplicity, we do not 
Secondly, we add contour lines of complex numbers of the same phase at 


integer multiples of 15 degrees (see Figure 4.1b). To that end, we change the 


intensity of the colour by multiplying it with a sawtooth function 

0.7 + (1.0 - 0.7) • {ip/iir/U) - L<p/(V12)J) ■ 

Because phase corresponds to hue, the points of such a contour line are all of 
the same colour. 

Finally, we add contour lines of complex numbers of the same modulus on a 
log-scale (see Figure 4.1c I. To that end, we change the intensity of the colour by 
multiplying it with a sawtooth function 

0.7 -h (1.0 - 0.7) • (log(r)/(7r/12) - [log(r)/(7r/12)J). 

Note that the contour lines of phase and modulus intersect each other ortho¬ 
gonally. The scaling factor l/(7r/12) in the sawtooth function for the modulus 
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Figure 4.1: Composition of domain colouring reference image: If we (a) represent 
phase by hue, |(b)| add contour lines of phase and |(c)| add contour lines of modulus, 
we obtain our domain colouring reference image, a variation of the enhanced 
phase portrait colour scheme of (Wegert, 2012[ Section 2.5). 


contour lines deliberately matches the scaling factor used in the sawtooth func¬ 
tion for the phase contour lines. Consequently, the regions enclosed by the 
contour lines of phase and modulus are squarish in appearance. 


4.2 Complex square root 


Recall the construction of a Riemann surface of the complex square root from 


Subsection 1.1 where we glued together its two branches at a branch cut along 


the negative real axis. 



(a) (b) 


Figure 4.2: Domain colouring of the two sheets of the complex square root 
The domain colouring of the two branches of the complex square root over a 
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square of side length 10 centred at the origin is shown in Figure 4.2 


On the sheet shown in|Figure 4.2a[ the complex square root takes values with 


negative real part (coloured green to blue). On the sheet shown in Figure 4.2b 


it takes values with positive real part (coloured purple to yellow). (The sheet 
shown in [Figure 4.2b| corresponds to the principal branch of the complex square 
root.) 

On both sheets, twelve contour lines of phase are visible, half as many as 
in the reference image. We can see that the phase of the complex square root 
function changes at half the angular velocity of its argument. 

Moreover, the discontinuity at the branch cut along the negative real axis is 
clearly visible. We also see that there is a smooth transition between the second 


(third) quadrant of Figure 4.2a and the third (second) quadrant of Figure 4.2b 


If we cut the two sheets along the negative real axis and glue the upper side 
of the cut of one sheet to the lower side of the cut of the other sheet, and vice 
versa, we obtain a Riemann surface of the complex square root. The resulting 
Riemann surface, produced with [Algorithm 2.1 


part as height function, is shown in Figure 4.3 
(multiview orthogonal). 


and [Algorithm 2.6 


(perspective) and 


using real 


Figure 4.4 



Figure 4.3: Domain-coloured Riemann surface of the complex square root in 
perspective projection 



(a) (b) (c) (d) 

Figure 4.4: Domain-coloured Riemann surface of the complex square root in 
orthogonal projection from left, front, right, and back (from left to right) 


Note that the self-intersection of the surface in [Figure 4T3| is only an artefact 
of using a height function to map the Riemann surface mesh from real four¬ 
dimensional to real three-dimensional space. Evidently, the two values of the 
complex square root at each point of the self-intersection do not agree: they are 
coloured differently, in green and purple, respectively. 
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In |Figure 4.4b we see the parabola that the real parts of the j/-values describe 
according to the equation — x = 0 when x takes values on the non-negative 
real axis. 


4.3 Folium of Descartes 



Figure 4.5: The folium of Descartes as a real plane algebraic curve 


The folium of Descartes is a classical plane algebraic curve of order three, 

C: f{x,y) = + y^-3xy = 0. (5) 

The cubic curve is nowadays called ‘folium’ after the leaf-shaped loop that it 
describes in the first quadrant of the real plane (see [Figure 4.5). It is named in 


honour of the French geometer Rene Descartes (1596-1650), who was among the 
first mathematicians to introduce coordinates into geometry. Originally, the curve 
was called fleur de jasmin since Descartes and some of his contemporaries, who 
were working out the principles of dealing with negative and infinite coordinates, 
initially wrongly believed that the leaf-shaped loop repeated itself in the other 



the y-value with the largest real part. 

We see that the first sheet carries y-values with negative real part (coloured 
green to blue). At the centre of the second sheet, we identify a zero of order 
two, which we can recognize from the fact that the colours of the colour wheel 
used in our reference image wind around it twice in the same order as in the 
reference image. It is the node of the leaf-shaped loop. The third sheet carries 
j/-values with positive real part (coloured purple to yellow). 
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(a) (b) (c) 

Figure 4.6: Domain colouring of three sheets of the folium of Descartes 


There are three branch cuts (discontinuities of hue) on the first sheet, six on 
the second sheet and three on the third sheets. We can see how the sheets of 
the Riemann surface are connected to each other along the branch cuts: First 
and second sheet are connected at the branch cuts of the first sheet. Second and 
third sheet are connected at the branch cuts of the third sheet. First and third 
sheet are not connected directly with each other. (Imagine how much harder it 
would be to read this off from Equation 5 ) 

Apart from the branch cuts, the map from x to y{x) is conformal (angle¬ 
preserving) on every sheet. We can see that the contour lines of phase and 
modulus intersect each other orthogonally on every sheet, as in our reference 
image. 

If we cut the sheets along the branch cuts and glue them together correctly, 
we obtain a Riemann surface for the folium of Descartes. 

The resulting Riemann surface, produced with|Algorithm !0]and|Algorithm 2.6 


using real part as height function, is shown in Figure 4.7 (perspective) and Fig¬ 


ure 4.8 (multiview orthogonal). 


Again, the self-intersections of the surface in [Figure 4.7| are only an artefact 
of using a height function to map the Riemann surface mesh from (real) four¬ 
dimensional to (real) three-dimensional space. 

[Figure TTj makes it obvious that cutting a Riemann surface into sheets by 
sorting y-values by real part may be the most straightforward but not necessarily 
the geometrically most appropriate method. Our Riemann surface of the folium 
of Descartes in large part appears to be composed of three copies of the complex 
plane (which looks like our reference image). Complications seem to arise only 
near the origin. 

If we look closely at [Figure 4.8b[ we may see how we obtain the real folium 
of Descartes (as a real plane algebraic curve) as the intersection of our Riemann 
surface mesh with the Re x-Re y-plane. The leaf-shaped loop is clearly visible as 
a hole in our visualization. One of the ‘complex planes of which the Riemann 
surface is composed’ is so thin that it is barely visible from this perspective. It 
is almost asymptotic to the ‘wings’ of the folium of Descartes (as a real plane 
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Figure 4.7: Domain-coloured Riemann surface of the folium of Descartes in 
perspective projection 


algebraic curve) in the second and fourth quadrant of the real xy-plane. 

Right below the centre of |Figure 4.8c[ we see two leaf-shaped loops in complex 
directions. Perhaps Descartes and his contemporaries were not entirely wrong 
after all to believe that the folium of Descartes has more than one leaf. Indeed, if 
we let x' = we discover that in the Re x'-Re y-plane the curve describes a 

leaf-shaped loop, which is exactly half as high as that in the Re x-Re y-plane (this 


also holds for the ‘wings’) and rotated into a different quadrant (see Figure 4.9). 


5 Conclusion 


We have discussed algorithms for the generation of a Riemann surface mesh of a 


plane algebraic curve (Algorithm 2.1) and its visualization as a domain-coloured 


surface (Algorithm 2.6) and their implementation using OpenGL and WebGL. 


The WebGL implementation combines floating point textures, multiple render 


targets, and a method due to Boubekeur and Schlick (2008) to replace the use of 
transform feedback and geometry shaders of the OpenGL implementation. While 
the generation of the surface takes noticeable time in both implementations, 
the visualization of a cached Riemann surface mesh is possible with interactive 
performance. This allows us to visually explore otherwise almost unimaginable 
mathematical objects. Sometimes the visualization makes properties of the plane 
algebraic curves immediately apparent that may not so easily be read off from its 
equation. It is possible to turn these domain-coloured Riemann surface meshes 
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Figure 4.8: Domain-coloured Riemann surface of the folium of Descartes in 
orthogonal projection from left, front, right, and back (from left to right) 



(a) (b) 

Figure 4.9: Leaf-shaped loops of the folium of Descartes in complex directions. 


into physical models using a full colour 3D printer. 
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